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Reaction paths and probabilities are inferred, in a usual Monte Carlo or Molecular Dynamic 
simulation, directly from the evolution of the positions of the particles. The process becomes time- 
consuming in many interesting cases in which the transition probabilities are small. A radically 
different approach consists of setting up a computation scheme where the object whose time evolution 
is simulated is the transition current itself. The relevant timescale for such a computation is the one 
needed for the transition probability rate to reach a stationary level, and this is usually substantially 
shorter than the passage time of an individual system. As an example, we show, in the context of 
the 'benchmark' case of 38 particles interacting via the Lennard- Jones potential ('LJ38' cluster), 
how this method may be used to explore the reactions that take place between different phases, 
recovering efficiently known results and uncovering new ones with small computational effort. 

I. INTRODUCTION 

The evolution of multi-particle systems arising in very diverse domains ranging from biology to material science is 
often governed by activated processes that occur with very low probability, but have a dramatic effect on the structure. 
Examples of physical phenomena controlled by such rare events are protein folding in biology, defect diffusion and 
crystal nucleation in condensed matter physics and cluster rearrangement in chemistry. 

Activated events are characterized by the fact that when they occur, their duration is rather short, e.g. of the 
order of the picoseconds in dense molecular systems, but that, in contrast, the typical time needed for them to start 
is very large. This separation of time-scale is the very definition of metastability; its origin may be energetic (high 
barriers) or, more likely in high dimensional systems, at least partly entropic (pathways hard to find). Even though 
physical times of activated events can be achieved in computer simulation using molecular dynamics, monitoring rare 
fluctuations is computationally expensive and still remains a challenge because of the long waiting time for the event 
to occur. To overcome this problem, several different strategies have been developed over the past years. A common 
idea behind these techniques is to bias the dynamics of the system in order to enhance the occurrence of reaction, 
or to use previous knowledge of the outcome of the reaction and to fix the endpoint of the trajectory. In order to 
implement this, in many cases one needs to know an order parameter able to discriminate between the reactant basin, 
the saddle regions and the product basin. 

A set of methods involving the direct sampling of path ensembles has been developed during the last decade^ 
The strategy is to restrict paths to the subset of reactive paths, those that interpolate between reactant and product 
basins. Examples of such methods are transition path sampling,—^ transition interface sampling.™ Another family 
of methods such as metadynamics^ and forward flux sampling 6 simulate the evolution in time of the system, and 
include some form of bias that guarantees that the reaction happens, but then of course the effect of the bias has to 
be discounted in order to retrieve the true statistics. Methods differ also in the way rate constants are estimated, 7 
and the ability to handle nonequilibrium steady states. 

Observation of rare events and exploration of energy landscapes of a mult i- particle system are indeed two intercon- 
nected problems. The requirement of knowing the target of the reaction, and in particular a relevant order parameter 
that parametrizes the path is problematic in physical situations where nothing a priori is known. Thus, a preliminary 
step prior to the calculation of finite-temperature rate constants using the sampling techniques mentioned above may 
consist in locating the saddle points of the energy and the corresponding energy minima using one of the eigenvec- 
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tor following methods described in the literature (e.g. the activation-relaxation technique, 8 Optim 9 or the dimer 
method 1 ^). The technique in all these methods relies on monitoring the eigenvalues of the Hessian matrix (the matrix 
of the second derivatives of the potential energy), in order to individuate stable or unstable directions, i.e. minima or 
saddles. A limitation of these methods is that energy saddles only correspond directly to the actual barriers for the 
dynamics if the temperature is very low, otherwise the entropic contribution to the dynamics becomes relevant. 

Another approach that may be used to explore the free-energy (i.e. finite-temperature) barriers between basin a 11 ' 12 
is inspired by Supersymmetric Quantum Mechanics, which has been used in the context of field theory to derive and 
generalize Morse Theory, precisely the analysis of the saddle points of a function. Transposing to statistical physics 
this formalism, one obtains a family of generalizations of Langevin dynamics, converging to barriers and reaction 
paths of different kinds, rather than to the equilibrium basins. The resulting method involves the evolution in phase 
space of a population of independent trajectories that are replicated or eliminated according to the value of their 
Lyapunov exponent The theory guarantees that by selecting trajectories having larger Lyapunov exponents in a 
specified way, the bias is just what is needed so that the population describes the evolution of the transition current, 
rather than that of the configurations, 14 as they would in an unbiased case. The advantage is that the convergence 
of the current distribution is much faster than the typical passage time. 

In this paper we shall present a more direct and elementary derivation of this dynamics without resorting to quantum 
theory (section [11]). Extending the more concise demonstration of Ref. I la. we will show how the modified dynamics 
precisely reproduces the evolution in the phase space of the probability current (or more precisely, the 'transition' 
probability current) between equilibrium basins, thus achieving a probability current sampling of the system dynamics. 
As the transition current evolves in time, it explores the different barriers, indicating which states are reached after 
different passage times. Starting from an initial equilibrium configuration, far from equilibrium phenomena are easily 
sampled, as the simulated current is an intrinsically out of equilibrium quantity^ 1 - In section IIII1 we discuss the actual 
population dynamic algorithm that may be used to simulate the probability current. Finally, in order to assess the 
performances of this probability current sampling, we present in section llVl an application to the structural transitions 
of a Lennard-Jones cluster, sampled under different physical conditions. We shall recover known results, and discuss 
some new ones, obtained in all cases at rather small computational costs. 



II. TRANSITION CURRENTS 
A. Overdamped Langevin dynamics 

Let us consider a 3D many-body system composed of N particles in contact with a thermal bath. We denote 
n, . . . , 7~3_/v the configurational degrees of freedom and V^({r"i}i = i...3jv) the interaction potential. In the limit of large 
friction 7 — > 00, where inertia can be neglected, the system evolves with a standard overdamped Langevin dynamics 
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where rji are independent Gaussian white noises of unit variance and rrij correspond to the particle masses. 

The probability to find the system at position r = (n, • • ■ r3jv) then evolves with the Fokker-Planck equation^ 

where we have introduced the Fokker-Planck operator Hpp- We can use the probability current Ji = 
~rm~ (j tB '^W' fr") P t° wr hc Eq. [5] as a continuity equation for the probability density: 

i 

For systems with separation of time scales, the dynamics can be split into two regimes. Starting from an arbitrary 
probability distribution, P(r; t) relaxes rapidly into a sum of contributions centered on the metastable states. At much 
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longer times, the rare transitions between the metastable states make P(r; t) relax to the equilibrium distribution. Two 
time scales can also be identified for the dynamics of the probability current. While the probability density rapidly 
relaxes into the metastable states, the probability current converges on the same time scale to the most probable 
transition paths between the metastable states. Then, the late time relaxation towards equilibrium corresponds to a 
progressive vanishing of the current, when forward and backward flux between each metastable state balanced Note 
that the same line of reasoning holds for non-equilibrium systems in which the forces do not derive from a global 
potential. In such systems, the probability current never vanishes and converges instead to its steady-state value. 

If one were able to simulate the evolution of the probability current, one would thus have all the knowledge relevant 
for the transitions between metastable states, while only having to simulate the system for relatively short time-scales 
(similar to the equilibration time within a metastable state). As mentioned in the introduction, simulating directly 
the transition current is the goal of this paper and we now derive a self-consistent evolution equation for Jj. 

Let us define the current operator Ji\ 

so that the probability current and the Fokker-Planck operator becomes 

J l = .hP (5a) 

H " = -Hw Si ' (5b) 

3 J 



The evolution of the probability current is then given by 

= JiP = -JiH FP P = - ]T J^JjP ■ (6) 
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where we have assumed that Hpp does not depend explicitly on time. Straightforward algebra shows that Ji-^r = 
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jf-Ji — fP Y and JiJj = JjJi which turns equation ([5]) into 



^-T.{lM-^3)p- (7) 

Using the expressions (|5a|) and (|5bj) for the currents and Fokker-Planck operator we obtain 

j% = —B.FpJi — - — - — jj . (8) 

3 J 

Note that the equations (|5a|) and ([7]) are not self-contained: the knowledge of P(r) is required to compute J^. On the 
contrary, ((5J depends exclusively on the current, and can readily be used to simulate Ji, without having to compute 
P(r) beforehand. The only condition is that the current distribution at the initial time Jf indeed derives from a 
probability distribution, i.e. is of the form: 



(9) 



This means that the initial current distribution should be such that the quantity Ai 

Ai = mi ei^T V J° (10) 

is a gradient, ^ L — j^r- A particularly simple initial condition is obtained if one assumes that P° is Gibbsean 

P° oc e k B TV in a region O, and zero elsewhere. Then, from Eq. ([9]), J® is zero everywhere except on the surface 
of fl, where it takes the form of a vector normal to the surface of J7, and with amplitude proportional to the Gibbs 
weight. 
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The evolution of current distribution given by Equation (|8|) , starting from an appropriate initial current J° converges 
to the stationary distribution of currents between metastable states on the same time scale as the usual Langevin 
equation converges to metastable-state. It is thus not necessary to wait for rare events to identify the transition path 
between the metastable states, an important improvement over standard MD methods. If there are several metastable 
states and transitions with different rates, the current distribution at longer times concentrates on the paths between 
regions that have not yet mutually equilibrated, and vanishes in transitions between states that have had the time to 
mutually equilibrate. 

B. Langevin dynamics with inertia 

In many physical situations inertia plays an important role and one cannot rely on overdamped Langevin equa- 
tions^ In such case, the stochastic dynamics of N interacting particles is given by 

n = Vi (11a) 
Vi = ~ TO ! rl ^~" ~JVi + \J 2m~ 1 ~fk B Ti]i , (lib) 

where («!,••• , v n ) = v correspond to the 3N velocities. 

In order to transpose the derivation of the current evolution given in subsection III Al to the inertial case, we have 
to consider the Kramers evolution equation P(r, v; t) = —HkP{y : v; t) instead of the Fokker-Planck equation^ 



dP 
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(12) 



As in the overdamped case, this can be written as a conservation equation where the probability current in phase 
space is given by 

J ri = Vi F(r, v; t) (13a) 

Once again, the current contains all the information about transitions between metastable states. There is however 
a conceptual difference: the presence of inertia makes it inherently difficult (and indeed, useless), to compute (| 13[) as 
it stands. The reason is that the phase-space current 
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is non-zero even in canonical equilibrium. For example, in a harmonic oscillator H = i [w 2 + r 2 ] , the phase-space 
current in equilibrium turns clockwise in circles around the origin. 

For this reason, the part of the current that corresponds to transitions between metastable states is screened by the 
large contributions of the currents within metastable states (see figure [1} . The probability current (|13a|) and (I13b[) 
does not really represent the transition paths between metastable states. 

We can however define a transition current 

Jl =J r+ _H__ 15a 

rrii ovi 

J* J kBT9P MSM 

J Vi = ^ ! ( 15b ) 

TOi OTi 

which has two interesting properties. Firstly, this current differs from the probability current by a divergenceless term 
and thus also satisfies the continuity equation P + V • J* =0. Fluxes out of a closed surface surrounding a metastable 
state are then the same for the probability and transition currents. The latter current thus contains the relevant 
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FIG. 1: Probability currents for an underdamped dynamics in a ID double-well potential. Left: Typical transition between the 
two wells. Center: Same trajectory in phase space. The probability current is dominated by oscillations in the wells. Right: 
Transition current J* in phase space. The transition path is sampled uniformly. Equilibrium contributions are not present. 



information about, for instance, transition rates. Secondly, the transition current vanishes in equilibrium, as can be 
checked by comparing (|14p and (|15|) . This current thus contains only the information relevant for the transitions 
between metastable states (see figure [lj and is not screened by the large 'equilibrium' currents within them. 

Using algebra similar to that of the overdamped case, one can showi£ (see the Appendix IVIj) that the reduced 
current evolves with 



dJ* 

~dt 



-HkJ 1 M J* 



with 



P 

p 

> J 1 I : 



(16) 



where the 6iV x 67V matrix M is given by 



M = 





1 d 2 V 



(17) 



Again, (|16l) is a self-consistent equation for the transition current and we shall now show how it can be simulated. 



III. USING POPULATION DYNAMICS TO SAMPLE THE TRANSITION CURRENTS 

The probability density P(r, v) is a scalar field that can be obtained by simulating many copies of the system, 
evolving their positions and velocities with the Langevin equation dill) , and constructing histograms. On the contrary, 
J' is a vector field and thus cannot be obtained in the same way: it also requires vectorial degrees of freedom. We first 
present a population dynamics that can be used to construct the evolution of the transition current in section UlI Al 
and then give the corresponding Diffusion Monte Carlo algorithm in section UlI Bl 



A. Stochastic dynamics 

For systems with many degrees of freedom, the direct resolution of the partial differential equation p^|) yielding the 
evolution of J* is not achievable numerically. In the same way as the Langevin dynamics (|11[) represents an alternative 
to the resolution of the Kramers equation, we can use a stochastic dynamics that generates the current evolution (|16[) 
numerically. 

The explicit construction of the dynamics was already presented in detail in a previous paperJ^ The idea is to 
proceed in two steps. Firstly, a Langevin dynamics (fTTj) is coupled to a carefully chosen evolution (see below) of a 
6iV-dimensional vector u = (u ri , . . . , u r3N , u Vl , • • • , u V3N ) in order to account for the vectorial degrees of freedom of 
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J*. Secondly, the transition current J' is given by the following average: 



6N 



J* = / f[ uF(r,v,u). (18) 

J i=l 

where f(r, v, u) is a joint distribution obtained from the Langevin dynamics, for instance by simulating a large 
number of copies of the systemJ^ 

In practice, the coupling between the vectorial and phase space degrees of freedom is obtained by the following 
population dynamics. One consider J\f copies of the system (called 'clones'), identified by positions and velocities r 
and v, which all carry a 6N dimensional unitary vector u. The dynamics of each clone is then as follows^ 

• r and v evolve with the standard Langevin dynamics with inertia 

• the vector u evolves with 

u = — M u + u(u'Mu) (19) 

• each clone has a birth-death rate a = — u^Mu. This is the only way the vector u influences the dynamics. 
The distribution of clones F(r, v, u) then evolves with 18 



i—l 



MijUj + ^(u^Mu) 



F - \JMuF (20) 



The first term of the r.h.s. comes from the Langevin dynamics (|11[) . the second one from the evolution of the 
vector (p~9j) and the last one from the birth-death events. The transition current is then given by flT8"|) . On can indeed 
check that taking the time derivative of the r.h.s of (fl~8| and using (|20|). one recovers the evolution of the transition 
current (fT6|) (see previous paper o 11 ! 12 for more details). 



B. Algorithm 

We first present the implementation of the population dynamics of the clones and then discuss the construction 
of the transition current. The m-th clone is denoted by its phase-space coordinates and vectors at time t, X™ = 
{r m (i), v" l (i), u m (t)}. We start with Af c clones whose positions and vectors are arbitrarily chosen. At every time 
step, the dynamics is as follows: 

1. All the vectors u m are rescaled to have a unitary norm. 

2. The positions and velocities of the clones are propagated using a leap-frog discretized Langevin dynamicsi 19 ' 21 

3. The vectors u m evolve with the (leap-frog discretized version of) the following dynamics: 

iiT = -M ijU f (21) 
Note that at after this step, the vectors are no more of unitary norm. 

4. For each clone one records w m = \\u. m (t + St)\\. 

5. We associate to each clone m a probability weight p m — M c w m / Wi and a random number e m chosen 
uniformly in [0, 1). The clone is then replaced by y m copies, with 

Vrn = YPm + f-ra\ (22) 

where \x\ is the integer part of x: if y m > 1, y m — 1 new copies of the m-th clone are made. If y m = 0, the 
clone is deleted and if y m = 1, nothing happens. The population size is thus increased by y m — 1 if y m > 1 or 
decreased by 1 if y m = 0. 
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6. After the step 5, the population is rescaled from its current size M^ = $^m=i Vm to its initial size M c , by 
uniformly pruning/replicating the clones. 

Steps 1 to 3 correspond to the propagation step of independent clones, whereas steps 4 to 6 correspond to selection 
steps. In particular, the steps 4 and 5 correctly represents the cloning rate a = -u'Mu of the previous subsection 
since ±\\u m (t)\\ = -u mt Mu"\ so that \\u m (t + St)\\ ~ exp(-<5iu mt Mu m ). 

The rescaling of the population at the step 6 can be done in many ways. For instance, one can pick a new clone 
at random M c times among the A/" c e clones obtained at the end of the step 5. We used an alternative approach 
that is less costly in terms of data manipulations: if 7V c e > M c , we kill W" c e — M c clones chosen uniformly at random 
among the Ml obtained at the end of the step 5. Conversely, if Ml < M c , we choose uniformly at random M c — Ml 
clones and duplicate them. Note that even though the clones evolve independently during the propagation step, their 
dynamics are correlated because the deleted clones are replaced by the duplicated ones at the selection steps. When 
the probability weights of the clones are all equal, we have p m — 1 and y m = 1. As a result, the population is left 
unchanged. Conversely, when the clone weights take distinct values, clones with small weights are likely to be replaced 
by those with large weights. 

There are many ways of implementing the resampling of the population (steps 5-6), well documented in the literature 
on Diffusion Monte Carlo . 22 i 23 In particular, it could be advantageous to do the resampling only every n time steps, 
where n is tuned to ensure ergodicity in phase space, i.e. to achieve enhanced sampling towards the unstable regions 
where saddles are located. 

Since the clones move in phase space with a Langevin dynamics, it can be surprising that they converge rapidly to 
the reaction paths, i.e. that they explore efficiently the transition states. This can however be understood by noting 
that their dynamics (without taking the averages (|18I0 is the so-called Lyapunov Weighted Dynamics— which is used 
to bias the Langevin dynamics in favor of chaotic trajectories. The clones will then tend to 'reproduce' favorably in 
the neighborhood of saddles, which are particularly chaotic regions of phase space, and to die in wells. This generates 
an 'evolutionary pressure' that helps the clone escape from metastable states and find the reaction paths. 

As mentioned before, this dynamics does not provide directly the transition current and one still has to construct 
the averages (fT8| . This can be difficult and clever methods to do so were discussed in the literature, for instance by 
Mossa and Clementi who studied the folding of chain of aminoacids. 24 The difficulty is connected to the well-known 
sign problem: large population of clones with arrows pointing in opposite directions cancel in the average but can 
numerically screen smaller asymmetric distribution that contains the information relevant for the transition current. 
One can however show that if one starts from a population of clones uniformly spread over a reaction path separating 
two metastable states and pointing in the same direction, the time taken for the sign problem to occur is of the order 
of the tunneling time through the barrier (see Appendix I VIII) . In the following we will always simulate much shorter 
times and omit the averaging steps to simply look at the distribution of clones. This distribution often suffices to 
locate the reaction paths, as can be seen on figured] for the double well potential. To extract further information, for 
instance regarding the reaction rates, we would need to do the averages, as was done in Ref. [2J, but this is beyond 
the scope of this article. 

IV. APPLICATIONS 

A. LJ38 cluster 

We now turn to the study of transitions between metastable states in the 38-atom Lennard- Jones cluster, a bench- 
mark model system that has been extensively investigated in the past . 19 ' 25 i 27 ~ — This system has a complex potential 
energy landscape organized around two main basins: a deep and narrow funnel contains the global energy minimum, a 
face-centered-cubic truncated octahedron configuration (FCC), while a separate, wider, funnel leads to a large number 
of incomplete Mackay icosahedral structures (ICO) of slightly higher energies. 

Although the lowest potential energy minimum corresponds to the FCC structure, the greater configurational 
entropy associated to the large number of local minima in the icosahedral funnel make this second configuration much 
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more stable at higher temperatures. As temperature increases, this system thus undergoes the finite-size counterpart 
of several phase transitions. First, a solid-solid transition occurs at T ss = 0.12^ when the octahedral FCC structure 
gives place to the icosahedral ones. At a slightly higher temperature, T s ; = 0.18-^-, the outer layer of the cluster 
melts, while the core remains of icosahedral structured This 'liquid-like' structure, also referred to as anti-Mackay 
in the literature, then completely melts around T s i — 0.35-^-d 

The numerical study of this system is challenging: global optimization algorithms have failed to find its global 
energy minimum for a long time 1 and direct Monte Carlo sampling fails to equilibrate the two funnels. The study of 
the equilibrium thermodynamics of this system required more elaborate algorithms such as parallel tempering^ - — 
basin-sampling techniques,— Wang-Landau approaches 32 or path-sampling methods . 19 ' 33 ' 34 

More recently, the dynamical transitions between the two basins has been studied following various approaches. The 
interconversion rates between the FCC and ICO structures have been computed using Discrete Path Sampling^ - — 
This elaborate algorithm relies on the localization of minima and saddles of the potential energy landscape, using 
eigenvector following, and then on graph transformation 3 ^ to compute the overall transition rate between two regions 
of phase space. To the best of our knowledge, this is the most successful approach as far as computing reaction rates in 
LJ38 is concerned. 38 However, the numerical methods involved are quite elaborate, require considerable expertise and 
have a number of drawbacks, all deriving from the fact that it is based on the harmonic superposition approximation 
and the theory of thermally activated processes. It thus requires any intermediate minima between the two basins 
to be equilibrated and this is only possible for small enough systems at low temperatures^ More importantly, when 
the difficulty in going from one basin to the other is due to entropic problems, as is the case for instance in hourglass 
shaped billiards, then the knowledge of minima and saddles of the potential energy landscape is not sufficient. 

Another attempt to study the transitions between the two funnels of LJ 38 relies on the use of transition path 
sampling. 33 Because of the number of metastable states separating the two main basins, the traditional shooting 
and shifting algorithm failed here, despite previous success for smaller LJ clusters^ 9 - The authors thus developed a 
two-ended approach which manages to successfully locate reaction paths between the two basins: they started from 
a straight trial trajectory linking the two minima, and obtained convergence towards trajectories of energies similar 
to those obtained in the Discrete Path Sampling approach. 33 Although the authors point out the lack of ergodicity in 
the sampling within their approach and the sensitivity on the 'discretization' of the trajectories, this is nevertheless 
a progress and the main drawback remains the high computational cost (the work needed 10 5 hours of cpu time) to 
obtain such converged trajectories. In contrast, the simulations we present below required less than 10 2 hours of cpu 
time. 



B. The LJ38 cluster and bond-orientational order parameters 

Before presenting our simulations results, we give some technical details on the LJ38 system and on the visualization 
of the different metastable states. The Lennard- Jones potential is given by the expression 

* ( I'" L : v) ><X 

j<k 

where r J = (r 3 x ,r^,r J z ) is the position of the j-th atom, rjk = |r J — r fe | is the distance between atoms j and k, e 
is the pair well depth and 2 1 / 6 <r is the equilibrium pair separation. In addition, all the particles are confined by a 
trapping potential that prevents evaporation of the clusters at finite temperature (i.e. particles going to infinity). If 
the distance between the position r of a particle and the center of the trap r c exceeds 2.25er, then the particle feels 
a potential |r — r c | 3 . LJ reduced units of length, energy and mass (cr = 1, e = l,m = 1) will be used in the sequel so 
that the time unit t = ayrnfe is set to 1. 

Rather than listing the 228 degrees of freedom of the atomic cluster, configurations are traditionally described using 
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FIG. 2: Short MD simulations were run to give an impression of the spread in the (Q6,Q4,E) space of each 'phases'. The 
simulation time was short enough that no tunneling between the phases was observed. The temperature was set to T = 0.15. 
The positions of the phases barely move in the (Qi,Q(,) plane when the temperature changes, although their spreading does. 
The kinetic energy however shifts when the temperature changes, and is roughly proportional to NkT where N is the number 
of degrees of freedom. 



the Qi bond-orientational order parameters**"! that allow to differentiate between various crystalline orders 

( 4 1 1 V /2 

Qi = 21TTN2 E fe'^)f > ( 24 ) 
V b m=-l ) 

where the Y/ m (#,</>) are spherical harmonics and 9jk,4>jk are the polar and azimuthal angles of a vector pointing 
from the cluster center of mass to the center of the (j,k) bond which connects one of the Nb pairs of atoms4* The 
parameter Q4 is often used to distinguish between the icosahedral and cubic structures, for which it has values around 
0.02 and 0.18 respectively^ Q 4 however does not distinguish between the icosahedral and the liquid- like phase and 
one thus often uses Qe, for which FCC, icosahedral and the liquid-like phase take values around 0.5, 0.13 and 0.05, 27 
respectively. To show the spread of the various basins in the (Q4,Qe) plane, we ran several molecular dynamic 
simulations, long enough to equilibrate within each basin but short enough so that one does not see tunneling (see 
figure [2|) ■ 

Although the whole temperature scale is interesting, the challenging part from a computational point of view is 
the low-temperature regime where ergodicity is difficult to achieve. Below, we show the results of our algorithm for 
three temperatures: T = 0.12e/feg, T — 0.15e//cs and T = 0A9e/kB that span the ranges around the solid-solid and 
partial melting transitions. 



C. Simulations 

Given the high dimensionality of the system, it is difficult to follow the evolution of all the coordinates of the clones 
in order to know if they have localized interesting structures. Instead, we proceed as follows: we plot the evolution of 
the average over the clone population of Q4, Qq and E as a function of the simulation time and we frequently store 
the positions and velocities of all the clones. 

If we see a plateau in Q4(t), Qe(t) and E(t), two cases are possible: either the clones have converged to a reaction 
path, or they are stuck in a metastable basin. In order to distinguish the two situations, we run an auxiliary short 
molecular dynamic simulation (without cloning) starting from the positions and velocities of the clones. The duration 
of this auxiliary simulation is long enough to observe relaxation into the metastable basins, but much shorter than 
the transition times. If the clones evolve away from the region they had populated in phase-space, we know they had 
found a reaction path and the auxiliary MD simulation converges to the metastable basins connected by this reaction 
path. If on the other hand the clones do not evolve away, we know that they had been stuck in some local basin. 
In such a case we can change two parameters to enhance the sampling of the phase space: the number of clones and 
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FIG. 3: Left and Center: Positions of the clones starting from the ICOm and FCC minima in the (Q4,,Q§) and (Qe,E) 
plans at T = 0.15. The clones starting from the icosahedral basin first find the barrier between ICOm and ICOam (black 
symbols, t = 1850). They then fall back in the ICO basin before finding a path that points towards the FCC funnel (blue 
symbols, t = 3500). Starting from FCC, the 600 clones find a path that leads toward the icosahedral funnel (green symbols, 
t — 1500). Right: We plot Qe as a function of time for the clones starting from the ICOm basins (red symbols) and the FCC 
basin (magenta symbols). Arrows indicates the time at which the snapshots shown in the left and center panel are taken. 



the friction 7 (see below for more details). The time step is always St — 0.01. Note that this procedure could be 
automated, but the way to do so is let for future work. 

In principle, any observable that can measure whether the population of clones splits in two separate sub-populations 
after a short Langevin dynamics would be suitable. If the clone population splits in two subpopulations with the 
same Q4, Qe and E, we may fail to detect the corresponding barrier. However, this coincidence would be extremely 
unlikely. 

Last, in addition to help us localize barriers, these short Langevin simulations allows us to explore the true dynamics 
close to a particular transition states. 



1. T = 0.15 

We first ran several simulations at T = 0.15, where the most stable state is the MacKay icosahedral minimum 
(ICOm) while the liquid-like phase (ICOam) and the FCC basin are metastable. 

Starting from the ICOm basin with Af = 200 clones and a low friction jSt — 10~ 3 , the clones rapidly find (t ~ 1500) 
a transition path to the liquid-like phase ICOam. Later on (t ~ 3500) an activated event bring the clones to another 
reaction path that points towards the FCC funnel. These times have to be compared with the transition time between 
the ICOm and FCC basins that was previously evaluated in the literature at roughly 10 7 .— Note that each barrier 
or path act as a metastable state for the cloned dynamics and it is by activation that the population jumps from one 
barrier to another. Running the same dynamics with a larger number of clones (Af = 600) tends to stabilize the first 
barrier so that one has to wait longer to see the transition to the second one. 

Starting from the FCC minimum with the same number of clones and at the same friction results in the clones 
rapidly going out of the FCC funnel and falling in the amorphous zone at the entrance of the icosahedral funnel.- 8 A 
reaction path is followed by the clones but not maintained. To stabilize this reaction path, we increased the number 
of clones and the friction. The effect of the former is mostly to slow down the dynamics while the latter allow the 
clones to populate the reaction path more uniformly. For Af = 600 clones and jSt = 1 , the population of clones indeed 
stabilizes the reaction path leading from the FCC basins to the entrance of the icosahedral funnel. The reason why 
we need more clones to stabilize this barrier than the ones in the icosahedral funnel is probably that the former is 
more flat and spread than the latter ones^l and thus requires a larger number of clones to be sampled uniformly. 

All these results are plotted on figure O in which we show the three basins ICOm, FCC and ICOam obtained from 
the initial MD simulations (see figure I2J and the positions of the clones in the {Q±,Qe) and (Qq,E) plans at different 
simulation times. 
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FIG. 4: Histograms made at the end of short MD simulations at T = 0.15 started from the clones positions at the times 
indicated by the green and blue arrows in figure [3] The gray-dotted regions correspond to the equilibrium MD simulations of 
the three basins ICOm, ICOam and FCC. Left: MD simulations started from the stationary structures found by the clones in 
the ICO funnel fall either back into the ICO basin or in a metastable basin around {Qi,Qa) = (0.3,0.05) that corresponds to 
an amorphous structure at the entrance of the ICO funnel. Center: The clones starting from the stable structure found in the 
FCC funnel fall either back in the FCC basin, or in a faulty FCC metastable state (blue rectangle) or in the ICO funnel, which. 
Both structures thus correspond to reaction paths. Right: Position of the surface atoms of a clone that has fallen in the faulty 
FCC configuration after the short MD. This figure was made using a Mathematica Spreadsheet that can be downloaded at 
http: //www- wales . ch. cam. ac .uk/~wales/make_f rames .nb 



To identify the various metastable basins connected by the clones, we ran several short MD simulations starting 
from the two long-lived plateaux (blue and green arrows in the right panel of figure [3]) . Histograms made from these 
MD runs are shown in figure [4] They show that the clones going out of the ICO basin find barriers toward the 
amorphous region at the entrance of the ICO funnel while the ones starting from the FCC minimum find a reaction 
path between the FCC basin and the icosahedral funnel. Interestingly, this path goes through a faulty FCC basin 
located around {Qa : Qq) = (0.12,0.45) that has been previously reported in the literature . 19 > 28 

The clones have thus found reaction paths pointing out of their starting funnels. The clones starting from the FCC 
basins find a reaction path that leads into the icosahedral funnel while the one started from the ICOm basin still 
remain in the icosahedral funnel. This could be explained by the fact that at this temperature ICOm is the stable 
state while FCC is only metastable so that the barrier ICO— ^FCC has to be harder to access than the one from the 
FCC side. Running short MD starting from the clones positions reveal intermediate metastable basins, either a faulty 
FCC or amorphous structures. 

2. T = 0.12 

This is the coexistence temperature between the ICOm and FCC minima. At such a low temperature, more and 
more secondary barriers play a role so that the transition between ICO and FCC becomes more and more complex. 
From the point of view of the time evolution of the transition current, this means that there are more and more 
metastable states for the clone dynamics. 

Starting from the ICOm basin with 600 clones and jSt = 10 -3 , we once again locate the barrier between ICOm 
and liquid-like phase ICOam. At this temperature this barrier is long-lived and we do not locate the one previously 
found at T — 0.15 that points toward the entrance of the icosahedral funnel. 

Starting from the FCC basin with jSt = 0.02, the 600 clones find several barriers that constitute a multi-step 
reaction path toward the icosahedral funnel. Once again, the larger the friction the longer the clones spend on 
intermediate barriers^ 

The position of the clones corresponding to the successive metastable barriers are shown in figure [5] Note that the 
typical times needed for locating the barriers are of the order of 10 3 , that is seven orders of magnitude smaller than 
the reaction times between ICOm and FCC, which is of the order of 10 10 i^ 
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FIG. 5: Left and Center: Positions of the clones starting from the ICOm and FCC minima in the (Qa, Qe) and (Qs,E) plans 
at T — 0.12. The clones starting from ICO find the barrier between ICOm and ICOam (black symbols, t = 8000). Starting 
from FCC, the clones find a succession of barriers that leads toward the icosahedral funnel (green symbols at t = 650, blue 
symbols at t = 1300 and cyan symbols at t — 2650). Right: We plot Qe as a function of time for the clones starting from the 
ICOm basins (red symbols) and the FCC basin (magenta symbols). 
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FIG. 6: The color codes correspond to MD simulations at T = 0.12 started on the green (left), blue (center) and black (right) 
arrows in figure [5] Left: Starting from the first stationary structure found in the FCC funnel, the clones relaxes mostly in 
the FCC basin and in the faulty FCC configuration shown in [4] Center The second barrier is close to the commitor between 
the ICO and the FCC funnel: the clone population relaxes almost equally in both funnel (57% of the clones fall back in the 
icosahedral funnel while 43% enter the fee funnel). Right Clones started from the barrier between ICOm and ICOam populate 
both basins. Note that the relaxation is much slower than for the other barrier because of the entropic nature of this barrier. 




Running short MD simulations starting from the clone positions on the barriers and constructing the corresponding 
histograms reveals various intermediate metastable basins in the ICO and FCC funnels (See figure [5]). The fact that 
Qi and Qq are not good reaction coordinates is confirmed in this figure: the first plateau (green points on figure [5]) 
seems to be after the faulty configuration when going from the FCC funnel to the icosahedral one but the MD starting 
from this barrier falls into the faulty configuration and the FCC basin, which seems to indicate that this barrier is a 
reaction state between the FCC and the faulty configuration. There is then a second barrier between the faulty FCC 
and the ICO funnel (blue dots in figure [5]). A last barrier leads to the amorphous region that separates the liquid-like 
phase and ICOm minimum. Note that in these regions the MD simulations are not very helpful because the transition 
between ICOm and ICOam has an entropic nature so that it is difficult to relax into the basins. The clones, however, 
successfully identify the barrier between these two states. 



3. T = 0.19 



This temperature is very close to the transition between ICO and liquid-like phase. As shown by free-energy studies, 
the barrier between the FCC and the ICO funnels is very low and the FCC basin is rather unstable.™ Clones starting 
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FIG. 7: Left: 100 Clones are started at T = 0.19 in the ICOm basin where they spend some time (first first time steps after 
t — 200, green dots) before locating the barrier toward the FCC funnel (first five time steps after t = 1700, blue crosses). 
Center: When the clones have found the barrier (t — 2000) a standard MD starts and relaxes as it should into the two funnels 
(black dots, five first snapshots after t = 2050). Right: Evolution of Qa(t). The cloning is stopped at t = 2000 and a normal 
MD follows. 

from the FCC basin do not stabilize on any structure because there is no proper 'rare barrier' and MD simulations 
starting from FCC immediately falls into the icosahedral funnel^ 

Starting 100 clones from the ICO basin at jSt = 0.01, they rapidly find a barrier connecting to the liquid-like phase. 
Later on, activated events lead the clones to locate a reaction path leading towards the FCC funnel. Starting MD 
from this barrier show that the clones relax into the FCC and ICO funnel. 

As mentioned above, it is hard for the clones to stabilize because the FCC funnel is barely metastable and the 
barrier crossed while going from FCC to ICO is rather flat at this temperature. It is thus quite surprising that they 
nevertheless manage to do so while starting from the ICO basin. If one starts from the FCC funnel, the clones almost 
immediately fall into the ICO funnel and then from there can locate the barrier again, but we were not able to stabilize 
the barrier when coming from the FCC basin. This might be due to the fact that clones stabilize reactions that take 
place on long time-scale (ICO— J>FCC), but not short-time relaxations (FCC— >ICO). 

D. T — 0.05: annealing the cloned system 
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Q 4 Q 4 

FIG. 8: Left: Result of a simulation run at T=0.12 with 600 clones, starting from the FCC minimum, with -ySt = 0.6 Right: 
Starting from the end point of the simulations at T=0.12, we run a standard cloning simulation at T = 0.05. After a time 
t — 1790 the 600 clones are still on the structure that had localized at T — 0.12, which is thus very stable. 

If one starts at such a low temperature from one of the various metastable basins, the clones remain trapped for 
a time longer than the simulation time. One can however use a temperature annealing to locate the barriers. If one 
starts the cloning simulation at T = 0.12 or higher, it is quite easy, as we saw above, to localize the barriers. The 
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temperature can then be decreased to T = 0.05 and the clones remains on the structure that were localized at a 
higher temperature (see figure [5]). 

V. CONCLUSIONS 

The algorithm we have discussed in this paper may be characterized as one that simulates the evolution in time 
of the current distribution, rather than that of the configuration. Because the time for the escape current to be 
established is often much smaller than the passage time itself, the method is able to find the transition paths very 
efficiently. 

The method has several attractive features: 

i) It does not require any previous knowledge of the relevant reaction coordinates. On the other hand, if an 
approximation of the reaction path is known a priori, one may always start the clones along this path, and they will 
populate the true current distribution in a shorter time. 

ii) Because the target of the dynamics is the reaction path distribution itself, one may perform simulated annealing 
in path space: first populating the reaction path corresponding to relatively high temperature, and then refining it 
to the lower, target temperature. Repeated annealing can also be used to locate several competing barriers in system 
with multiple reaction mechanisms. 

Hi) The reaction current vanishes between mutually thermalized regions^ This is why at longer times, the system 
converges to the barriers that take longer to cross, irrespective of whether they are of entropic or energetic nature. 
This may be an advantage in cases in which the energy landscape is not in itself dominant, but rather the multiplicity 
of paths dominates. 

iv) The construction of the transition current and the cloning algorithm also applies for non-equilibrium systems 
where the forces derive locally, but not globally from a potential, such as a system with leads at the edges having a 
potential difference. Reaction paths between non-equilibrium metastable states, which cannot be described in term 
of a free energy, may be studied in the same way. The only difference is that the average (|18p does not vanish in the 
long time limit and converges instead to the steady-state transition current. 

Note that since the method does not require the knowledge of the reaction coordinate, it could be used efficiently 
in systems with competing reactions where one does not know a priori the end points of the reaction paths. This 
would for instance be particularly interesting when studying the crystallization of suspensions of oppositely charged 
colloidsi^JS 

In principle, the reaction time may be expressed directly in terms of the (unnormalizcd) reaction current. It may 
also be recovered from the weights carried by the clones, which may possibly be achieved from importance sampling 
in a Lyapunov-weighted ensemble of trajectories^ However, the method, as it stands does not allow one to calculate 
the reaction time with great precision, due to the exponential nature of the timescale. Further work is required in 
this direction. 
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VI. APPENDIX A 

We report here the derivation of the time evolution equation (|16p for the transition current in the underdamped 
(i.e., Kramers) case of section IPTBl 

The classical Kramers probability current J, presented in equations (|15[) , can indeed be written, as in the Fokker- 
Planck case (Eq. (|5ajl). in the operatorial form J(r, v, t) = JP(r, v, t), where the components of the current Kramers 
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operator are 



9 fo-i 7 d 1 dV 



m.j 9rj 



(25a) 
(25b) 



In equations (| X5[) of section III Bt a transition current J* has been introduced, which can be expressed in turn with 
operators as 

J* = J f P = (j + t) P (26) 

where J is the Kramers current operator reported above, giving the usual phase-space current J, and the 'transition' 
operator 

1 d 



TV 
TV 



(3mi dvi 
d_ 
dr.; 



(27) 



TP is a divergenceless term. As already remarked in section IIIB1 the transition current still satisfies the continuity 
equation 



^- = -V, lV • J 
at 



(28) 



thanks to the divergenceless of TP. We have introduced here the phase-space divergence V rjV = (V r , V v ) 

As in section Til Al we proceed now in deriving explicitly the time evolution equation of J*. Multiplying both sides 
of the continuity equation above (indeed identical to the Kramers equation (|12l) ) by the transition current operator 
leads to 



- dP 

J*^r = -J*Vr,v • J*P 

at 



(29) 



On the l.h.s. the transition current operator can be commuted with the time derivative. The r.h.s. of (|29l) can be 
rewritten with commutators as 



J*V r v • J* = V r . v • J ( J + T 



J, V r v J 



T, V r v • J 



(30) 



using (|2"9"]l and the zero divergence property of T. 

Resorting to definitions of the current operator J and the transition operator T given in (|25a|) . (|25bl) and (j27|) . 
explicit expressions for commutators in (|30[) can be recovered with straightforward algebra: the term J, V r>v • J 
gives 



J f . \/ r v " J 



p = -J2 s ™ [ J -» 



7 
/3m, 



d Vi P 



1 d 2 V 

rm dridr a 



v,P 



(31a) 
(31b) 



while 



T, V r ,v * J 



can be expressed as 



P = J2(/3m i )- 1 5 ia {8 ri + 7 d Vi )P 

i 



(32a) 
(32b) 
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Inserting now (|31"j) and (|32j) in the r.h.s. of Eq. (|29|) yields 

(4 + 2V )(V r , v • J)P = (V r>v • J)J^ - ^ 5 io (J Ka - 7 (/3m l )- 1 a„ l P) 

+ ^(/3m l )- 1 ( 5 la (5 n + 7 9 Ul )P (33a) 

i 

(X a + r«.)(V r ,v • J)-P = (V r ,v ' J) J* + y^iafrJo. " TT^P + — ^-J U ) 

i <9 2 y 

a» 4 p (33b) 



/3m! 2 dridr a 



that can be recasted as 



leading to equation (jI6 



VII. APPENDIX B 



In our study of the LJ38 system, we carried out simulations of the stochastic dynamics of the clones but did not 
explicitly make the averages (|18j) that would yield the transition current. We argued in section Ull Bl that if one only 
looks for the reaction paths, then these averages are not necessary. In this appendix, we illustrate this claim on a 
simple Id example that also allows us to discuss several aspects of the clone dynamics, namely the metastability, the 
finitcness of the clone population and the role of the initial condition. 

We thus consider a system undergoing an underdamped Langevin dynamics 



x = v; v — —71; — m 1 V'(x) + y/2'ykT '/mr] (35) 

where V(x) is a potential with two barriers, plotted in figure [9] We ran the clone dynamics for 2000 clones starting 
in the left well and carried out the averages (|18p. 




FIG. 9: Left: Plot of the potential V(x) = x(-39 + 240x + 15x 2 - 138a; 3 + 20x 5 )/120. Right: The green crosses are the position 
of the 2000 clones after a time t — 400. The black arrows correspond to the averages (|38[) and indeed point tangentially to the 
reaction path. The color code and contour lines corresponds to the value of the Hamiltonian H(x,v) = v 2 /2 + V(x) 
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To do so, we constructed an approximate density from the positions and vectors (x l , v l , it*, ) of each of the 
J\f c = 2000 clones: 

1 M, 

Fnum{%, v,u x ,u v ) = j^/^2S(x - x l )S(v - v l )8{u x -u l x )S(u v -u l v ) (36) 

c i— 1 

In principle, the S should be Dirac functions but for practical purposes we replaced the one acting on the phase space 
coordinate by the bell-shaped function 

S n (x, v) — — exp [ 2 2 | if x 2 + v 2 < w 2 and S n (x,v) — otherwise (37) 

where w — 0.1 and Z is a normalization constant. Finally, using (|36|) and (|37|) in ()18j) . we construct the transition 
current from the numerical data by computing 

1 ^ 

3 T (x,v) = — rVu'^^-^tJ-Uj) (38) 

on a grid every = = 0.15 and plot the resulting vector if its norm is larger than 10 -3 . For visualization purposes, 
we plot in the figures the vectors 5 times longer than they really are. 

We started a simulation with 2000 clones in the left well, around x = mi ~ —1.9, with unitary vectors (u x ,u y ) 
pointing at random. The temperature is set to kT = 0.09 and the friction to 7 = 1.5 so that the mean first passage 
time across the barrier isH 



V7 2 /4+|V"(M 1 )|- 7 /2V V "M 

where Mi is the first maximum of the potential Mi ~ — 1. The results of the simulation after a time t = 400 are 
plotted in figure M The clones have already populated the barrier4£ As can be seen, the averages (|3"5|) along the 
reaction path are non-zero and result in vectors tangent to the reaction path, pointing toward the left well. 

At later times, two processes take place, roughly on the same time scale. Firstly, more and more clones come 
back from the central well to the left one. Their vectors u are always tangent to their trajectories, but can be 
pointing toward the left or the central well with equal probability. Indeed, if (p(t), q(t), u(i)) is a possible trajectory 
of the system, then so is (p(t),q(t), — u(t)). As a result, the averages (13"8"|) may cancel out at large times, when the 
subpopulations of clones whose vectors u point toward the central and the left wells balance. This is how numerically 
the transition current is supposed to vanish at large time (another possibility being that all the clones leave a region 
of phase-space, because of finite population-size effects). 

Secondly, some clones reach the barrier leading to the right well and duplicate, which results in populating the 
second reaction path. Since the clones did not have time to fall in the right well and cross back the barrier towards 
the central well with vectors u pointing in the opposite direction, the average (|38|) does not cancel along this reaction 
path. 

Both effects can be seen in figure [TU] the clones populate both barriers; the average ((38)) cancels out over the first 
barrier but not yet over the second one. This shows that the clone dynamics do locate the barriers and remain on the 
reaction paths even though the transition current may average out when the two wells separating a barrier equilibrate. 
This thus validates our approach to locate the reaction paths in the LJ38 system. 

Note that if one wishes to study quantitatively the transition current, two modifications would need to be done to 
our algorithm. First, the initial condition should not be taken at random but constructed as proposed in section Hi Al 
Second, rather than simulating all the clones in parallel while maintaining their population constant, it may be 
advantageous to run them sequentially, starting one run for every offspring of every clones, as is done for instance 
with the 'Go with the winner' methods4£ The constrain on the total population being fixed indeed affects the 
metastability of the clone dynamics and increases finite size effects. For instance, if there are J\fi and A/2 clones on 
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FIG. 10: Left: t = 997 the clones populate both barriers. The arrows average out along the reaction path between the left and 
central wells which have equilibrated, whereas the transition current is still present between the central and right wells Right: 
at t — 1590 the clones only populate the reaction path between the central and right wells. Since the simulation had enough 
time to equilibrate the involved wells, the average (fl8|) cancels out. 



the same reaction path, with vectors pointing in opposite directions, both populations grow exponentially with the 
same rate. Now, if the total population is kept constant, then the smallest sub-population disappears on average 
exponentially fast. Using a 'Go with the winner' method would palliate this drawback but would result in additional 
computational costs difficult to estimate beforehand. 
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